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Abstract. The unique fluctuation-dissipation theorem for equilibrium stands in 
contrast with the wide variety of nonequilibrium linear response formulae. Their most 
traditional approach is "analytic" , which, in the absence of detailed balance, introduces 
the logarithm of the stationary probability density as observable. The theory of 
dynamical systems offers an alternative with a formula that continues to work when the 
stationary distribution is not smooth. We show that this method works equally well for 
stochastic dynamics, and we illustrate it with a numerical example for the perturbation 
of circadian cycles. A second "probabilistic" approach starts from dynamical ensembles 
and expands the probability weights on path space. This line suggests new physical 
questions, as we meet the frenetic contribution to linear response, and the relevance of 
the change in dynamical activity in the relaxation to a (new) nonequilibrium condition. 
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1. Introduction 

An open system in contact with a large environment is in stationary equilibrium for a 
given reduced scale of description when at each moment its condition realizes minimal 
free energy. The total entropy production is then zero and the evolution is time-reversal 
invariant. Small disturbances break the stationary behavior and after some time an 
eventually new equilibrium is established mainly via dissipative effects. The comparison 
with the original equilibrium is the domain of linear response: since about 60 years 
now PQ[2], f° r classical, quantum, open, and closed systems basically the same formula 
relates the response to a small perturbation with an equilibrium correlation function. 
This formula (later called Kubo formula) exactly picks up the physical interpretation in 
terms of dissipation, hence the name fluctuation-dissipation theorem [3]. 

Response theory continues making sense for stimuli to nonequilibrium systems, 
where entropy is being produced already before the perturbation. Over the years 
various types of linear response formulae have indeed been obtained for nonequilibria, 
including rather diverse but specific models such as in climatology and for glassy or 
coarsening dynamics. For better orientation the present paper identifies some common 
traits between the various approaches to put them in a more unified framework. Suppose 
for example that one wants a numerical code for predicting the response of an out-of- 
equilibrium systems, i.e., without actually perturbing the system. Without a clear 
picture of the available theoretical results, the choice of the method, if any, would be 
limited and possibly suboptimal. To answer requests like these, we aim at classifying 
the features of several "extended fluctuation-dissipation theorems" for nonequilibrium, 
highlight briefly their strengths and weaknesses, and the eventual relations between 
them. 

Reviews and classifications of nonequilibrium response already exist. The paper 
by Marini Bettolo Marconi, Puglisi, Rondoni, and Vulpiani has reviewed response in 
statistical physics in the light of recent fluctuation relations [4]. Seifert and Speck have 
introduced a classification of some fluctuation-dissipation theorems into three classes [5] . 
Chetrite and Gupta present a more mathematical view [6]. The present work includes 
the case of deterministic dynamics and attempts a concise classification of linear response 
for nonequilibrium. We will conclude that there are three main classes. The first class 
of formulae can be derived from the Kubo-Agarwal formula, which itself starts from a 
Dyson-expansion of the perturbed Markov semi-group. A second approach originates in 
the theory of dynamical systems and can be applied when the stationary distribution 
is not smooth; it gives rise to a numerical algorithm which in fact also applies for 
certain stochastic dynamics, as we show. The third class is much more probabilistic 
and treats noise as an important observable in the linear response. The synthesis is 
there provided by introducing the excess in dynamical activity in a second "frenetic" 
contribution to the traditional Kubo fluctuation-dissipation formula. In that sense, 
the name fluctuation- dissipation theorem (even "extended") is not fully suitable for 
nonequilibrium systems, as e.g. their return to stationary nonequilibrium is not uniquely 
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characterized by dissipation [TJ. On the mathematical side, our classification is to work 
either with generators (for Markov evolutions) or to work with weights on path space 
(mostly still limited to classical dynamical ensembles). Referring to those, we call it the 
analytic approach (Section [2]) versus the probabilistic approach (Section [3]). Remarks 
complementing those in previous sections, and conclusions, are found in Section HI 

2. Analytic approach 

The framework we consider is that of Markov dynamics (no memory) on regions of 
M. d . That includes deterministic dynamical systems, being essentially first order in 
time. It also includes jump processes as described with Master equations, but we use 
here the language of diffusion because it reduces naturally to that of deterministic 
dynamical systems in the limit of zero noise, realized here with a state-independent 
diffusion constant D — > 0. For even much greater simplicity we choose overdamped 
diffusions where the velocity field 

x t = v(x t ) is given by (1) 

v(x t ) = F(x t ) + V2D£ t 

for standard white noise £t [S] . In Section 13.21 we also treat (underdamped or inertial) 
Langevin processes. In all these cases it makes sense to speak about the so called 
backward generator L, working on observables (Heisenberg picture). The expectation 
of an observable A at time t > is then given by 

(A(x t )) = J p(dx) (e tL A) (x) (2) 

when at time zero the states are distributed with probability p, possibly singular with 
respect to the reference volume element dx on state space. We can also abbreviate 
that as d(A(t))/dt = (LA(t)). When no confusion can arise, we continue to write 
A(t) = A(xt) for the (most often random) value of the observable at time t > 0, and 
(B(0)A(t)) for the time-correlation between observable B at time zero and observable 
A at time t. 

For the overdamped diffusion (CQ) the generator is 

(LA) (x) = (F ■ VA) (x) + DAA (x) (3) 

A perturbation changes the drift F — > F h = F + hF\ where, for simplicity, we avoid 
inserting a time-dependence in the small amplitude h, applied at all times t > 0. This 
leads to a perturbed backward generator L h = L + hLi with L\A = F\ ■ VA. The change 
in expectations at times t with respect to what we had for the unperturbed dynamics 
follows from (j2J): 

(A(t)} h - (A(t)) = J p(dx) (e tLh - e tL )A (x) (4) 
To linear order in h, 

e tLh - e tL = [ e sL {L h - L) e (t ~ s)L ds + 0{h 2 ) (5) 
Jo 
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yielding a linear response 
d(A(t)) h 



dh 



(Mt)) h - (A(t)) 
h 

= X (t) = [ dsR(t,s) 



= lim 

h=Q h^O 



including a susceptibility x{t) as the integration of a response function of the form 



R{t, s) = I p{dx) {e sL U e {t ~ s)L A) (x) 



(6) 



(7) 



More generally, when applying a time-dependent perturbation h s at time s > we also 
have 

mt)) h 

5h s h=o 

We mostly restrict ourselves to the case where p is stationary: the response depends 
only on the time difference r = t — s and we can write 



R(t,s) 



s <t 



(8) 



R(-, 



p(dx) L l e rL A (x) 



(9) 



as the central object of study for the linear response of stationary Markov evolutions 
within an analytic approaclijj. 

The direct reading of the right-hand side of (Q, further discussed under Section 
12.21 is that e rL A evolves the observable A for a time r and then L\ acts on the result 
to evaluate it in state x. However, we start in the next section with the more frequent 
approach of acting on p. 

2.1. Acting on probabilities 

In this section we focus on manipulations with the stationary probability distribution p. 
The basic step from (Q is partial integration, which means that it is assumed here that 
p has a smooth density with respect to the reference volume element, p(dx) = p(x) dx. 
In many cases that appears to be a reasonable physical assumption when the level of 
description is mesoscopic to macroscopic, independent of whether the system is driven 
or not (see [9]). 



2.1.1. Kubo-Agarwal formula Assuming a smooth density p we have that © can be 
rewritten as 



R(t) 



dx p{x) 



P 



x) e rL A (x) 



(10) 



where the adjoint C\ is defined by J dx(Cip){x) A(x) = f dxp(x) (LiA)(x). Adjoints 
are forward generators of the time-evolution on densities, as appears e.g. in Master 

| We can of course also use ((9|) for discrete jump processes. For each pair x, y of states, the perturbation 
enters as a modification of jump rates w(x,y) — > w[x,y) + hw±(x,y). The matrix L\ to be put in (j9|) 
has elements w\{x,y) \ix^y, and W\{x,x) = — *}2yWi{x,y). 
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equations. For the diffusion ([3]) the adjoint of L is the operator of the Fokker-Planck 
equation (d t p = Cp, see [8]) 

Cp = -V • (Fp) + A(Dp) (11) 
and C\p = — V ■ (F\p) so that (flOil takes the form 

i?(r) = -([V • *i(0) + F 1 (0) ■ V logp(O)] A{t)) (12) 
which is a specific realization of 



mr) = (^(0)A(r) 
P 



(B(0)A(t)) (13) 



with observable B Note that in general the stationary expectation (5) = 

because J dx C±p(x) = from the normalization of p. Applications of that Agarwal 
formula [10] in practice meet the difficulty of needing to know the density p (which is 
usually not available) and the details of the dynamics for C±. It is thus a result (from 
partial integration) on a formal level. 

Formula (fl~3l) is associated to equilibrium, see formula (2.5) in [TT]. Of course we 
have only used that p is smooth. It is easy to verify that in the case of detailed balance 
with conservative forces F = —VU at temperature D — 1/(3 we obtain the Kubo formula 
[2] for the linear response relation in equilibrium. Indeed, say for ([3]) with perturbation 
Fx = — VV around equilibrium p oc exp[—/3U], we get 

-It = -pw -VU + AV = (3LV (14) 
P 



C p (0) A(t)\ = P (LV(0) A(r)> = ±{V(0) A{r)) (15) 



so that 

It 

\ p I dr 

which is the (classical) Kubo formula. The last equality has used that under detailed 
balance (LV(O)A(t)} = (V(O)LA(t)) (see further details in Sec. I2"X3"|) . In a way, the 
Agarwal formula (|T3|) repeats Kubo's original derivation while stopping short before 
specifying p. 

Also others have re- found the Agarwal formula, such as in Theorem 2, formula? 
(2.22) and (2.23), of Chapter 2 in [12J. Weidlich gives a quantum version: his equation 
(2.17) replaces L\p |[p, Hx] with the commutator of the perturbing Hamiltonian 
Hx [13] • Hanggi and Thomas in [T3] find the Agarwal formula in their equation (3.12) 
for time-dependent processes. In the review [1] formula (2.70)-(2.73) is the Agarwal 
formula. We also find these formulae such as f JT2|) involving log p presented in the book 
by Risken, formula (7.10)-(7.13) in [8], and as formula (7) in [15J. There also re-started 
the emphasis on logp as "generalized" potential. 

2.1.2. Information potential The formula that Falcioni, Isola and Vulpiani [16J derived 
for tiny displacements of the initial condition along a unit vector e, 



R(r) = (e-Vlogp(0)A(r)) (16) 
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(see it also as formula (3.13) in the review [4]) corresponds to the case of diffusion 
with constant perturbation F\ = — e (impulsive constant force in the direction — e) in 
( Tl2i) . Again, as there is no explicit dependence on the noise level the formula can be 
readily tried for chaotic dynamics with smooth stationary density [17J. The function 
X = — logp is sometimes called the information potential. In fact, that potential gets 
a prominent place in various works on nonequilibrium linear response that follow the 
Hatano-Sasa formalism [18] , such as in the more recent [19] or [20] . 

Other formula? focusing on X can be found in the works by Prost, Joanny and 
Parrondo [20] and by Speck and Seifert [5]. Consider the stationary density p h for 
the perturbed dynamics with generator L h and the linear response p h (x) — p(x) = 
hpx(x) + 0{h 2 ). From stationarity C h p h = we get C\p + Lp\ = 0. For the Agarwal 
formula ( TTBl we need 

B = ^ = (17) 



On the other hand, 



d h l h = d h \ogp h = \im 



p h - p pi 



h-+o hp h p 

(dh is understood in h = 0) so that for inserting in ( fl3l) . B = C{pdhl h ) / p. We conclude 
R(r) = [ dxC(pd h l h )(x)e rL A(x) 



dx p(x) d h X h (x) Le rL A(x) 

— [ dxp(x)d h X h (x)e TL A(x) (19) 
dr J 



or 



R(r) = ^-(d h l h (0)A(r)) 



(20) 



This formula with the special emphasis on the presence of the information potential 
appears in [5] (but with a different time-derivative), where X is called stochastic entropy. 
The result of [20] starts from expanding the Hatano-Sasa identity, which effectively 
makes a special choice for the observable A. There one imagines the process and hence 
also the X to depend on a family of parameters A and one asks how the expected value 
of A = d\X at A = A* changes under a small change in these parameters around a given 
A*. The answer is provided by (|20~1) . which is equation (5) in [20J. 

As already noticed for (Tl3|) a disadvantage of incorporating the information 
potential into the correlation function of ( !T6|) and (|20|) is that it is generally unknown 
and not directly measurable. Not only do we not know the stationary density, but also 
no clear physically relevant interpretation in terms of heat, work or thermodynamic 
potential has been found for X. On the plus side, a clear advantage of ( 1T6|) and 
( 120]) is that no detailed information on the dynamics is needed and that one can use 
parametrized forms of the stationary density p to put in X. Typically, (quasi-) Gaussian 
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approximations are tried for p to make formulae (I12p . (1161) . (1201) more explicit and 
practically useful; see e.g. [211 EI2J [221 122]. Note that there is no reference in formula 
( 1121) to the noise strength except through the stationary p. If one replaces the p in the 
B of (1131) by a Gaussian, its variance will effectively reflect the noise level but needs to 
be fitted. 

Finally, recent works on time-dependent processes with feedback rewrite the 
positivity of the entropy production in terms of the expected information potential 
(relative entropies) which appears useful for understanding work relations, see e.g. 
[211 [25]. At any rate, ( 120]) is the climax of the analytic approach for smooth probability 
densities, preserving the Kubo form ( f!5l) most faithfully. 



2.1.3. State velocity There is a relation between the information potential X and the 
state-space velocity u. For diffusion processes with generator ([3]) this state velocity is 

u = J -^ = F - DV\ogp = F + DVT (21) 
P 

as can be readily checked from the expression for the stationary probability current 
j p appearing in the Smoluchowski equation d t pt + Vj^ = expressing conservation 
of probability p t - Of course this probability velocity needs not be related to physical 
currents, but an interesting observation writes a nonequilibrium response formula as a 
co-moving (equilibrium) fluctuation-dissipation theorem. 

In a stationary process two quantities A and B have time-translationally invariant 
correlations (B(0)A(t)) = (B(—t)A(0)). As we mentioned above, the Kubo formula 
(or, the fluctuation-dissipation theorem) holds under detailed balance, i.e., for an 
unperturbed evolution which gives rise not only to a stationary but also to a time- 
symmetric distribution of trajectories. Nonequilibrium means breaking time-reversal 
invariance. We introduce the operator L* that generates the time-reversed motion to 
describe for example 

(B(-t)A(0)) = fp(dx)A(x) (e tL *B) (x), t>0 (22) 

One can alternatively write J p(dx) (Lf)(x) g(x) = J p(dx) f(x) (L*g)(x), from which 
we see that L*g = C(pg)/p. For the overdamped diffusion processes that we have 
considered so far, from Q3J) and PJ, L* = L - 2u • V = -FV + DA + 2DV logpV. 

Note now that the Kubo formula ( lisp would follow from the Agarwal formula when 
in (USD we take B = -(3L*V. Indeed, 

(L*V(0)A(t)} = (V(O)LA(r)} = ±(V(0)A(t)) (23) 

Detailed balance corresponds to L = L* (i.e., (B(—t)A(0)) = (B(t)A(0))) and 
corrections to the Kubo formula will thus arise from choosing in ( 1131) 



B(x) = ±_ v (x) = -L*V(x) ^— V (x) (24) 



for some function V. In other words the antisymmetric part L 2 L V in the formula §M 
for B will be responsible for violating the Kubo formula and will show the nonequilibrium 
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aspect. Moreover, L — L* = 2u ■ V relates to the state velocity (12~T|) . Substituting these 
relations via (j2l|) in f lT3|) . we find that for diffusions fl3]) (with D = 



The equation (1261) shows that the equilibrium Kubo form gets "restored" when 
describing the system in the Lagrangian frame moving with drift velocity u. The 
passage to the Lagrangian frame of local velocity d/ds — > d/ds — u ■ V "removes" 
the non-conservative forcing from the formula, as explained in [151 [26l 127] . Still, if we 
do not know the stationary p, the formula ( 1261) implies a statistical average over the 
unknown vector u. 

2.2. Acting on the observable 

We go back to the original (Q and we move to an algorithm focusing on the evolution of 
the observable A rather than that of the density p. Various contributions to the theory of 
dynamical systems start exactly from that formula for the formulation of linear response 
as there still no assumption is needed on the smoothness of the probability p. 

The main player now is L\e tL A in (JTj), and L need not to commute with L\. 
Typically the generator of the perturbation is Li — F\ • V for states x G M. d . We 
are thus interested in obtaining a useful relation for Ve tL A. In other words, we need to 
start from 



Ruelle treated a formula with this form for deterministic dynamical systems [281 ITT] - 
recently applied especially in studies of simple models for climate response [221 EHl EQ] . 
Indeed (1271) is especially suited in case the stationary probability law is strange in the 
sense that its lack of smoothness forbids further partial integration to go back to the 
Agarwal formula (Tl3|) . In that same context however, it is useful to further split (TSTJ into 
two parts [TT] , one describing fluctuations along the stable directions of motion and one 
parallel to unstable directions (those with positive Lyapunov exponents). Informally you 
would expect that the vector Ve rL A (x) can be given, at least for large r, in its natural 
components exp[rAe(x)] e-VA(x) for local Lyapunov exponents corresponding to the 
various (stable and unstable) directions e. Such a decomposition is natural for stationary 
distributions p which belong to the class of Sinai-Ruelle-Bowen (SRB) measures; these 
have a density in the unstable (expanding) direction (positive Lyapunov exponents) so 
that there a further partial integration remains possible towards an Agarwal formula 




(25) 



abbreviated as 



R(t, s) = p±(V(s)A(t)) - 0([u(s) ■ W(s)]A(t)) 



(26) 




(27) 
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(TT3"j) . That is exactly what is suggested in the hybrid form of Section 2.3 of [22J. 
See [3TI 132] for an introduction on SRB-measures with their physical interpretation. 

2.2.1. Algorithm for deterministic and stochastic systems The numerical evaluation of 
linear response for chaotic dynamics has recently been studied by various groups. In 
particular, Abramov and Majda have developed new computational approaches based 
on the theory of SRB-measures [29]. We simplify here the presentation to introduce 
an algorithm that works also for stochastic systems. To exploit ( 1271) . no knowledge 
of p is required, but we need to get a hold on V x e rL A (x) = V x (A(x T )) X0=x where we 
have emphasized that the differentiation is on an evolved quantity with respect to the 
initial condition x. To solve the problem that V (state derivative) and L (generating 
the time-evolution e tL A(x)) need not to commute, we unfold the formula a little further. 
A practical numerical tool for estimating V x A(x T ) in deterministic dynamics (where x T 
is uniquely determined from x = x) was already presented in [33J, but the following 
numerical method is more efficient for steady states. 

Consider first an evolution in discrete time n = 0, 1, . . . (but with a parameter e 
that will allow a continuous time limit e — > 0), 

x n+ i = g n (x n ) = x n + v n (x n ), xq = x (28) 

with v n (x) = e[F(x) + £ n ] 

for £ n anything stationary in time, including possible "noise" depending on the time 
n, but that does not depend on the state x. This (£„) is considered frozen so that for 
its given realization we take (1281) as a deterministic dynamics. The main point is that 
V ' x {A(x n )) x = (V x A(x n )) x where the (•) averages over the "noise" (£ n ). We then need 
to deal with V x A(x n ) where x n depends on the initial state x through ( 1281) . By the 
chain rule, applied recursively, 

V x A(x n ) = (VA)(x n ) ■ V x {g n -i o...og ) (x) 

= (VA)(x n ) ■ G n _i (x) 

= (VA)(x n ) (Vp n _i)(:r n _i) ■ V x (g n - 2 o...og ) (x) 

= (VA)(x n ) (V<7 n _i)(ar„_i) • G n _ 2 (x) = . . . (29) 

where V A is a 1 X d row vector (easily computed at all times) and Gk(x) = V(gk°gk-i ° 
. . . o g ) (x) is a d x d matrix obeying the recursive relations 

G n (x) = Vg n {x n ) ■ Q n -i{x) (30a) 

Go(ar) = Vg (x) (30b) 

Note now that for each time k, the analytical form of the matrix Vg^ = I + eV-F does 
not depend on k or on the "noise" (£ n ). Hence, 

V x A(x n ) = (VA)(x n ) ■ [I + £ (VF)(vi)l 

• [I + e(VF)(x n _ 2 )] •••[! + e(VF)(x)\ (31) 

The continuous time version Gt{%) can be obtained in a suitable limit n — > oo for G n (x) 
with time t = ne and time step e = t/n. In this limit (1281) returns the evolution 
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&t = F(x t ) for the deterministic case, while for stochastic equations an additional 
suitable rescaling to reproduce e.g. Gaussian noise with £ t is needed. At any rate, 
from I + e(VF)(i) — e e ^ F ^ x ' we get formally 



V x A{x T ) = (VA)(x T ) -Texp 



ds(VF)(x 8 



(32) 



where each x s depends deterministically on x for frozen (£ r , < r < s), and T indicates 
a time-ordered integral. The formula (1321) or its discretization (1311) is ready to be inserted 
into ( 1271) where the (•) will first average there over the possible "noise" and then, with p, 
over the initial condition x. That provides the main algorithm for the analytic approach 
working on the observable instead of on the probability distribution. 

A similar expression can be found in the Appendix B of Ref. [33] for the case of 
deterministic dynamics. In that work, however, the estimate of (VA)(x n )&t{x) was 
performed with an adjoint scheme, by integrating numerically backward in time the 
final value (VA)(x n ) with the equation 

d t (VA)(x t ) + (VA)(x t ) ■ (Vg)(x t ) = (33) 

(In [33] they wrote the equation for the column vector rather than for the row vector 
VA, and a transpose of (Vg) was used). Although that scheme is equivalent to our 
( l3TT) . it requires a CPU time 0(n 2 ) for a perturbation active during all n iterations, 
as opposed to the 0{n) matrix multiplications ( 12.2. ID for estimating matrices for 
k < n. Our scheme is exploiting stationarity: the propagator from time n — k to time 
n is the same as Gk from time to k. In transient regimes we would lose this property 
and we would also need 0{n 2 ) operations. 

A further study of such a numerical algorithm is contained in [22] [29]; in particular 
Appendix A of [29] explains the derivation above (without the generalization to 
stochastic evolutions). For questions of dealing in a similar context with the impact 
of stochastic perturbations, we refer to |34j . 



2.2.2. Numerical illustration To illustrate the numerical scheme ( )29|) for estimating 
Ruelle's linear response formula (1271) in a simple context, we consider a set of equations 
introduced in biology to describe circadian cycles, that is the periodicity of biorhythms, 
for Drosophila [351I36J. The state space has d = 5 dimensions, with states x = (Po, P±, 
P 2 , P/v, M). The dynamics couples the concentration M of mRNA with those of four 
types of proteins, written as Po, Pi, P 2 , and Pn in Ref. [36], where one can find the 
details of these equations. Denoting 



P 
~K 



P 



P + K 

we have equations of motion x 
dP 



F(x) of the form 



dt 
dl\ 
dt 



k?M 











'Pi 






+ 










P\ 








V2 






"3 



p\_ 



+ ^4 



(34) 

(35 a) 
(356) 




Figure 1. Evolution of the concentration of mRNA and of the four proteins of the 
model from random initial conditions, according to (|2.2.2|l . 



dP, 



dP 2 




'Pi' 




"P 2 ~ 




"P 2 " 


dt 


= "3 


[K 3 _ 


- I/ 4 


[K 4 _ 


K d 



hp 2 + fc 2 pj 



JV 



iV 



dt 

dM 



hP 2 - k 2 P N 



N 



- V n 



M 
K 



(35 c) 
(35 d) 
(35 e) 



3, u 3 = 6, u 4 = 3) nMh" 1 



1.5, K 4 = 2) nM 



with parameters as in previous papers, 

(u s = 0.5, v m = 0.3, v x = 6, u 2 

(K m = 0.2, Kj = 2,K 1 = 1.5, K 2 = 2, K 3 

(k s = 2, h = 2, k 2 = 1) h" 1 , n = 4 

Integration of (12. 2. 2D was performed with a simple Verlet scheme with a discrete time 
step e = 0.0025. Starting from random concentrations, the model reaches quickly a 
cyclic regime, as shown in Fig. [TJ 

We check the response in the mRNA concentration M to a change of the rate 
ks ks(l + h), starting from the steady state, i.e., a random phase of the cycle. 
The dynamical perturbation introduced by that change is HF\ = (h M ks, 0, 0, 0, 0), 
since a change in ks affects only the evolution equation of variable Pq, see (135 ap . The 
observable A(x) = M has a gradient (VA)(x) = (0,0,0,0,1), which is coupled to 
the perturbation via the "propagator" matrix G n . With the setup of Sec. 12.2.11 we can 
estimate G n by a sequence of matrix multiplications. Given an evolution in discrete time 
g(x) = x + F(x)e, we are interested in matrix Vg(x) = I + VF(x)e, or in coordinates 



(Vy)i 



5ij + dFi/dxj e. The rows of the matrix (VP) arc 



'iV 



V2 



Po 



'IV 

K 2 



0, 0, k s 



(36 a) 



Pi 
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Figure 2. Linear response function R(t) of the mRNA density M to an impulsive 
change in the parameter kg — > fcg(l + ho) as a function of time and response 
6(M(t)) /Sho computed with ho — 10~ 3 . Units are nM vs. hours. 




(0,0,^,-^,0) 
0,0,0,^ ' 



pr, 



M 



K 



(366) 

(36 c) 
(36 d) 
(36 e) 



m J M 



where derivatives with respect to P are denoted as [j^] 



P - - (p+Kf and similarl y 

for 9^" • This matrix with x Q = (P , Pi, P 2 , Pn, M) yields CVg)(x ) = I+(VF)(x )e, 

coinciding with G . Iteratively, Gi = (Vg)(x 1 ) ■ G , and so on. 

By sampling many trajectories, in parallel for perturbed and unperturbed dynamics 
starting from an initial condition from the steady state, the estimate of 5M(T)/5h has 
been obtained by setting a small constant h for r > and by evaluating the time 
derivative of ((M(r)) ft - (M(r)))/h (we used h = lO' 3 ). The response function instead 
has been calculated with (|2"9j) and (}2"7]) . and well overlaps with the response, as shown 
in Fig. [2j Asymptotic stability in this model (36] probably favors a good performance 
of formula (|27|) . We computed both the perturbed dynamics and the unperturbed 
fluctuation formula just to show graphically that the results are the same, but of course 
there is no additional understanding of the process. Obviously, computing the response 
only with unperturbed simulations would be convenient as long as the convergence of 
the method is good. We postpone a more detailed study of the efficiency of this and 
other numerical schemes to a future work. 
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3. Probabilistic approach 

Linear response makes sense in general only within a statistical theory. That is to say, 
sensitive dependence on initial and boundary conditions can create strong and relevant 
effects beyond linear order on the microscopic scale while macroscopic linearity remains 
valid. For estimating the mobility we do not investigate the microscopic particle's 
individual motion at specific times and how it changes under an external field, but 
we ask for a spatio-temporal averaged current and that includes noise. In fact, for 
stronger microscopic chaoticity we expect the statistical approach to be more relevant. 
While various physical variables can show chaotic behavior in their time-evolution, their 
spatial or temporal averages will typically have a much smoother behavior, see also [9] 
and chapter 6.2 in [37] answering the so called van Kampen objection. 

3.1. As noise gets important 

In coarsening dynamics of low temperature spin systems, or in spin glasses, a very long 
transient regime may exist towards equilibrium, i.e., the nonequilibrium is not imposed 
by external gradients and the dynamic equations are actually equilibrium ones. It 
is in this context that an extensive response-literature has been produced during the 
last decades. We mention briefly some results from 2003, whose focus was mostly to 
develop zero-field algorithms, in other words, exploiting fluctuation-response relations 
to estimate numerically responses from fluctuations in unperturbed dynamics. Works 
of Chatelain [38], Ricci-Tersenghi [39], and Crisanti and Ritort [10] introduced new 
schemes for computing the response of the system from correlations in the spins. They 
were followed by Diezeman [H] and by Lippiello, Corberi, and Zannetti [12} H31 HI] . The 
derivations of these results cannot be reduced in our brief scheme and various forms of 
response relations have been obtained. However, the one by Lippiello et al., a specific 
case of ( l4Tj) below, emerged as having clear physical significance beyond numerical 
usefulness, because it matched the form of an earlier study by Cugliandolo, Kurchan and 
Parisi [15] , who proposed a response relation for autocorrelations in Langevin equations. 
With the present notation (pQ), perturbation h s x and observable A(x) = x, that would 
read 



showing an "asymmetry" term in addition to the form of the Kubo formula. The 
derivation of (I3"7j) is based on the equality (x t £t) = 2Df3R(t, s), hence it is centered on 
the presence of noise £ t . This was new, certainly with respect to derivations presented 
in previous sections where eventually noise was just an aspect of the dynamics, not 
fundamental for the derivations. An equation similar to ([57)1 appears in a more recent 
study for nonequilibrium steady states |46j, based on a path-space formulation by 
Harada and Sasa [I7J, see their Appendix B. The Harada-Sasa approach has really 
pioneered the path-space approach of the next subsection. The consistency between the 



R{t, s) 



2D 



1 




) - {(F{x t )x s ) - (F{x 8 )x t )) 



(37) 
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approaches of Refs. [ISJ H21 HZ] indeed shows that these belong to a general framework 
with new significant physical content. In the following section we discuss the path- 
space formulation embracing these results and we mention the physical interpretation 
that goes beyond merely dissipative aspects. Noise becomes visible then as dynamical 
activity, ruling the time-symmetric fluctuations. 

3.2. Path space approach 

The origin of dynamical ensembles is the projection of a microscopic Hamiltonian 
dynamics on the dynamics of reduced variables. That Mori-Zwanzig projection [48| [49] 
originates from making a physical partition of the phase space, depending on the physical 
situation at hand, in which each microstate X is mapped (many-to-one) to a reduced 
state x(X). That induces noise in the reduced dynamics, for example on the mesoscopic 
level of description. It is then natural to consider dynamical ensembles, i.e., probability 
distributions on path space where a path refers to a trajectory on the mesolevel. Such 
was already the approach of Onsager and Machlup starting dynamical fluctuation theory 
in [50]. These trajectories, under certain limiting conditions (e.g. via a weak coupling 
limit or via adiabatic elimination), can be described via first-order equations, in which 
case we meet the Markov processes of the previous Section [2j Dynamical ensembles can 
however also describe non-Markovian processes describing important memory effects, 
see e.g. [5T] for the application of a response relation in a visco-elastic medium. 

Paths are trajectories uj = (x s ,0 < s < t) in state space, say looking at the 
states x s in the time-interval [0, £]. For the sake of simplicity we characterize the 
unperturbed (perturbed) process by the probability weight P(to) [P h (u))] with respect 
to some reference du. The mathematical idea is to turn perturbed expectations 
{A(t)) h into unperturbed ones (A(i)P h (u>) / 'P(u)\. Defining the relative action U{u) = 
\ogP(u)/ P h (u) and splitting U{u) = [T(u>) —S(u)]/2 into a time-symmetric T(w) and 
a time-antisymmetric S(u>), we get 



where Ti and Si are the linear contributions in h of T and S, respectively, around h = 
(where U = 0). This formula is general but the point is that in physical systems with 
local detailed balance jl?] we have a good physical understanding of the path functions 
T and S E31 EH ES]- Let us first look under global detailed balance. We already see 
here that since the observable A(t) — A(Q) is time-antisymmetric (ignoring momenta), 
under detailed balance 





(386) 




{A(t)) h - (A(t)) = {A{t) - A(0)) 



h 




(39) 
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because ([A(t) - A(0)])7i(w)) = and (A(t) «Si(w)) = -(A(O)Si(w)) by time-reversal 
symmetry of the reference equilibrium process. The result f )39|) should equal the Kubo 
formula ( fl5l) . and indeed it does as will become clear. When the perturbation is from a 
potential —h s V(x s ), it is well established that the time- antisymmetric S(u>) = Si(u>) 
must be the path-dependent entropy flux into the environment as caused by the 
perturbation potential V [561 EZ]- If the environment is at uniform temperature, for 
small constant h for times t > the change in entropy is S(u) = /3h[V(t) — V(0)], namely 
dissipated energy h[V(t) — V(0)] divided by temperature. The first term on the right- 
hand side of (138 cp has thus the same form one finds in equilibrium (1391) , but in general the 
unperturbed process is in steady nonequilibrium and S(u) = (3h[V(t) — V(0)] is an excess 
of entropy production with respect to that already generated by the nonequilibrium 
dynamics. Hence, the derivative 8(A(t)S(u)) /5h s = (3-^(V(s)A(t)) also mimics the 
equilibrium Kubo formula. 

We now turn to the second term on the right-hand side of (138 cp . Eq. (138 cp must 
be true also for a constant A, for which the response is zero. We thus get 



If -^-T\{u) equals a state function B(s), then, since ^-(V(s)) = (LV(s)), one deduces 
that B(s) = (3LV(s), arriving at 



For diffusion processes (JTJ with generator and temperature D = 1/(3 we see that 
indeed LV(x) = (F- W)(x) + DAV(x) is a state function. The same is true for Markov 
jump processes (for specific derivations of (HTT) using stochastic calculus, see [HI]). 
The function LV quantifies the time-symmetric volatility of V under the unperturbed 
dynamics, also named "frenesy". It is related to the "dynamical activity" in discrete 
systems, where it quantifies the change in escape rates; T{oS) is the excess in dynamical 
activity over [0, t\. We refer to [53HM1I55] for more details and examples. Thus, there is a 
new contribution next to the first term in (I4ip (taking one half of the Kubo formula ([151) 
and referring as usual to dissipation). In particular, its form now clearly deviates from 
response formulae such as (1201) . One then wonders in what sense the decomposition in 
formulae like ( T4T1) is intrinsically natural. We think it is, as that splitting makes the first 
term time-antisymmetric and the second term symmetric in time s, which corresponds 
to the dissipative and the reactive part of the susceptibility, respectively. 

Note that in (1411) one computes averages without needing an explicit knowledge 
on the stationary probability. On the other hand some knowledge is required on the 
dynamics sitting in the generator L. It is still not clear how much kinetic information 
is truly needed and how practical that gets. We summarize the general physical idea in 
the formula 




(40) 



R{t,s) 



2 cb 



(V(s)A(t))-'-(LV(s)A(t)) 



(41) 




(42) 
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where EntrP'^w) is the excess in entropy flux over the time period [0,t] and Esct '^ (co) 
is the excess in dynamical activity due to the perturbation. The physical challenge is to 
learn to guess or to find that Esc[ ' T '(a;) from partial information on the dynamics. For 
the overdamped dynamics of e.g. (I4ip it means to consider the expected rate of change 
of the perturbing potential V under the original dynamics. The response relation out of 
equilibrium is no longer a fluctuation-dissipation relation but a fluctuation-dissipation- 
activation relation. In fact, the formula can now be turned around and from measuring 
violations of the fluctuation-dissipation relation one obtains information about the 
active forces [5Tj . 



3.3. Inertia! case 

We open a separate subsection on the inertial case of Langevin dynamics because its 
linear response is much less discussed in the literature despite the obvious interest for 
example for models of heat conduction. The main point however is that the ideas 
summarized under (13. 2p and ( 142|) remain unchanged. 

For states (q,p) = (q 1 , q 2 , . . . , q n ] p\p 2 , • • • ,p n ) € M. 2n of positions and momenta 
we attach standard white noise Q to each 1 < i < n, with constant strength D l and a 
friction coefficient 7* to model heat baths at temperature D % /Y = T l : 

q l =p l 

f = F\q) - 7 V + h t ^r + (43) 

The forces F l can contain a nonconservative part but are confining when we want a 
stationary regime where the particles typically reside in a bounded region. We already 
inserted the perturbation V(q) with small time- dependent amplitude h s for s > 0. The 
linear response is given by formula (17) in [52] : 



-Zw<{(^' )F ' Mm 



d I dV , , A / , 



+ Y.{-^M.u£m)} (44) 



3 

The first sum again corresponds to the dissipative part from the entropy fluxes in the 
reservoirs at temperatures T\ The remaining sums give the frenetic contribution. As 
not recognized yet, formula (jUJ) can still be rewritten in a similar way as done in ( l4Tj) . 
Supposing D % = j/fl we must replace there V — > V/j = p ■ V q V/j, 



R(t, s) = £±<p. ■ V g V(s)A(t)} - ^(L(p s ■ V q V)(s) A(t)) 



(45) 
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in which the generator L for (j43p now reads 

Lf(q,p) =p-V q f + (F- 7 p) • V p / + DA p f (46) 

and is the underdamped version of ([3]). Note however that the dissipative last term 
DApf does not contribute in ( 143]) and that these formula? do not work in case D l = 
for an i with dV/dq l ^ 0. In that case the best alternative is probably to apply the 
algorithm of Section 12.21 

An important application of f|44|) is to the modification of the Sutherland-Einstein 
formula relating transport coefficients such as mobility with fluctuation quantities such 
as the diffusion constant [581. 



4. Further remarks and conclusions 



4-1. Nonlinear responses 

There is also a growing number of works on higher-order terms around equilibrium. In 
fact part of the book by Evans and Morriss is devoted to that [59]. Other references 
include [60J, Section 10 in [57] or the more recent (H EU [62], [63]. One typical start 
is the fluctuation symmetry in the distribution of the entropy flux, transient from the 
reference equilibrium system as also explained in [61] for thermostated dynamics. 

As was emphasized in [65J the main point is probably not to be able to write formal 
expansions and formulae, but to find useful structures and unifying interpretations. 
We will not deal with that here, except for mentioning one particular general relation 
between second order and first order terms that has been largely unnoticed [66J, and an 
instance of which has appeared as identity (25) in [65] . 



2 



dh 2 



(A(t)) 



= P^-([A(t)-A(0)][V(t)-V(0)]) h 



(47) 

h=0 



The A equals A up to flipping the sign of the momenta. All terms are explicitly 
expressed as correlation functions in the equilibrium reference process (indicated with 
the superscript (-) ). The derivation of (l4"7j) uses linear response around nonequilibrium 
for non-state functions (i.e., the correlation function in the right-hand side). 



4-2. N on- state functions 

In case of observables that depend explicitly on time or are functions of the trajectory 
over several times (such as products A 1 (t 1 )^4 2 (^2) • • • A n (t n )), some of the formula? must 
change. The basic techniques remain however in place. A typical application is how 
the heat depends on a change in parameters, e.g., for estimating nonequilibrium heat 
capacities [67J. Heat is not a state function but varies with the trajectory and the 
applied protocol. 
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4-3. Note on effective temperature 

A traditional approach to violation of the equilibrium fluctuation-dissipation theorem 
is to imagine an effective temperature in the otherwise unchanged Kubo formula [68J. 
That idea has had most success with mean field type systems, but it has remained 
more unclear how the effective temperature can provide a consistent and general tool for 
realistic systems. Still today it serves as a paradigm for interpreting experimental results, 
see e.g. [69]. One possible approach for the future would be to associate an effective 
kinetic temperature to the ratio between the frenetic and the entropic contribution in 
(142|) . see Appendix A in [53]. 

4-4- Transient case 

Lots of attention have been devoted to the linear response behavior for relaxational 
processes. We mentioned some of that in Section [3TTJ, but there is no way to be complete. 
For example, the interest in ageing and glassy dynamics has much stimulated the search 
of modified fluctuation-dissipation relations [70]. Here we emphasize that the analytic 
approach and in particular the methods of Section I2T1 loose their simple structure when 
the unperturbed reference is time-dependent (and not stationary as was assumed). On 
the other hand, formulae like ( l4Tj) are unchanged when one is not starting at time zero 
from the stationary distribution, but one needs to take the average in the unperturbed 
transient regime. Also causality is automatically verified there as the second term in 
( )4T|) equals the first term when s > t. This unification in expressions for the transient 
and the steady regime is only natural as even the stationary regime is physically and 
ultimately but a very long transient. 

4-5. General properties of linear response 

Some formal properties of linear response go basically unchanged from equilibrium to 
nonequilibrium contexts. For example, sum rules and Kramers-Kronig relations only 
depend on the presence of an underlying Hamiltonian dynamics for the total system 
plus environment or on causality. Of course, specific expressions can differ but there 
is physically nothing new, see also [7_lJ [72]. Nevertheless, there do exist essential 
differences. We have already alluded to the fact that the name "fluctuation-dissipation 
theorem" is no longer so appropriate because of the importance and complimentary 
character of changes in the dynamical activity (and not only in the dissipation). In 
fact, the word "fluctuation" also becomes less correct as the time-correlations in the 
expressions for R(t, s) no longer express a symmetrized time-correlation. To make that 
point, let us evaluate the equilibrium Kubo formula ( fl5|) ) when the perturbing potential 
V equals the observable A: 

(V(t)) h - (V) = -(V(0)[V(t) - V(0)]) = \{[V{t) - V(0)f) > (48) 

which is the variance (or the "diffusion") of the displacement V(x t ) — V(xq). That is 
sometimes called a generalized Einstein relation. In fact all linear transport coefficients 
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for equilibrium can be expressed like that, as was understood already in 1960 by 
Helfand [73]. For nonequilibrium linear response, that relation including its positivity 
gets violated, see also [581 E]- I n the wor ds of the previous subsection, negative 
effective temperatures become possible in nonequilibrium. The origin lies in the frenetic 
contribution, e.g. the second term on the right-hand side of ( 14 1 p . which can overrun the 
first dissipative term. 

4-6. Outlook 

Response theory is primarily about predicting the reaction of a system in terms of 
its unperturbed behavior. As we mostly have in mind response in time, that involves 
temporal correlations. Therefore, dynamical fluctuation theory can be expected to be 
most prominent. (See [71] for the distinction between static and dynamic fluctuation 
theory: statics looks at deviations around the law of large numbers at single times, for 
example for the average over many copies of the system, while dynamical fluctuations are 
around the law of large times, for deviations around time- averages.) Around equilibrium, 
dynamical fluctuations are governed, just as static fluctuations, by the entropy and 
dissipation functions. That further enables connecting different response coefficients, 
such as in Onsager reciprocity or via Maxwell relations. Nonequilibrium makes a more 
drastic difference between static and dynamical fluctuations. In particular, no useful 
connections between different types of responses have been discovered for nonequilibrium 
processes. We think of the analogue of relations between compressibilities, heat 
capacities and conductivities. In our opinion, that challenge will require finding 
experimental access to quantities like the dynamical activity which is complementing 
entropic characterizations in the description of nonequilibrium processes. 

4-7. Conclusions 

We have presented a concise guide (with plenty of references) to the multi-faceted 
world of linear response for systems out of equilibrium. From there we have discovered 
similarities but also some missing pieces that otherwise would not be evident. One 
such piece was the extension to stochastic systems of Ruelle's formulation, which we 
have introduced, together with an efficient algorithm. That also allows for numerical 
calculations where an explicit knowledge of the density of states is not required, as 
opposed to other (analytic) formulations that we have described. However, the resulting 
linear response formula contains a correlation function whose physical meaning is not 
very clear. 

In contrast, a rich physical picture emerges from a probabilistic approach based 
on path space weights, where the stationary distribution is also not needed. Besides 
the quite different mathematical apparatus compared to the more standard analytical 
approach, the probabilistic way indeed emerges as the one that currently offers 
more relations with dynamical fluctuation theory: one has to study also how the 
system correlates with the activity of the perturbing potential, a time-symmetric 
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quantity complementary to the time-antisymmetric fluxes of entropy. The combination 
fluctuation-dissipation therefore does not suiffice to characterize the linear response of 
nonequilibrium systems. 

Acknowledgments M.B. thanks M. Colangeli, A. Vulpiani and J. Wouters for useful 
discussions. 
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